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PREFACE 


As  part  of  its  Project  RAND  research  program,  RAND 
engages  in  basic  support  studies  in  mathematics.  The 
present  Memorandum  is  concerned  with  computer  techniques 
combining  several  methods  of  approximation. 


SUMMARY 


Dynamic  programming,  successive  approximations, 
extrapolation,  and  smoothing  are  used  in  this  RAND 
Memorandum  to  treat  ill-conditioned  systems.  Numerical 
examples  are  given. 
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DYNAMIC  PROGRAMMING  AND  ILL-CONDITIONED  LINEAR  SYSTEMS 

1.  INTRODUCTION 

A  system  of  linear  equations 

(1.1)  Ax  =  b 

is  said  to  be  ill-conditioned  if  A  is  almost  singular. 

The  computational  solution  will,  in  general,  require 
some  additional  devices  since  the  numerically  large  elements 
in  A-1  amplify  small  changes  in  b,  and  small  numerical  errors 
incurred  in  the  course  of  the  calculation,  to  the  point 
that  they  overwhelm  the  significant  elements. 

In  a  number  of  situations,  particularly  in  those 
of  engineering  and  physical  origin,  we  have  some  a  priori 
information  concerning  the  solution  x.  For  example,  we 
may  know  that  |  |x  —  c|  |  is  small,  where  c  is  a  known 
vector,  or  we  may  know  that  the  components  of  x  are 
monotone  increasing, 

(1.2)  £  Xg  £  •••  1  Xjj. 

In  the  first  case,  we  can  use  this  additional 
information  by  considering  the  new  problem  of  minimizing 
the  quadratic  form 

(1.3)  (Ax  —  b,Ax  —  b)  +  X(x  —  c,x  —  c)  , 

where  (x,y)  denotes  the  usual  inner  product  and  X  is 
a  parameter  to  be  chosen  adroitly.  This  technique  has 
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been  used  by  Twomey  til  and  Phillips  [2] . 

In  the  second  case,  we  could  consider  the  problem 
of  minimizing  the  quadratic  form 

(1.4)  (Ax  —  b.  Ax  —  b)  +  X(y2x,y2x), 

2 

where  v  x  represents  the  vector  whose  i— th  component 
is  xi  —  2xi_i  +  xi_2*  with  xQ  **  x^  =  0. 

In  this  paper  we  shall  confine  our  attention  to  a 
consideration  of  the  problem  posed  in  (1.3)*  There  are 
several  features  of  novelty  to  our  treatment.  In  the 
first  place,  we  shall  employ  a  dynamic  programming 
algorithm,  cf.  (3!,  1 43  ,  to  determine  the  minimum  of 
the  expression  in  (1.3)*  Secondly,  we  shall  combine  this 
with  successive  approximations  In  x  and  extrapolation 
In  X.  Examples  will  be  given  of  the  efficacy  of  these 
methods. 

2.  DYNAMIC  PROGRAMMING  APPROACH 

To  apply  dynamic  programming  to  the  minimization  of 
the  quadratic  form 

(2.1)  RjjU)  0  (Ax  -  b, Ax  -  b)  +  x(x  -  c,x  -  c), 

we  observe  that  the  minimum  value  is  a  function  of  b, 
and  indeed  a  quadratic  function  of  b. 

Consider  the  more  general  problem  of  minimizing 
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(2.2) 


Vx> 


X^xl  "  cl)4 


+  X(xg  -  c2)2  + 


Mxj*  “ 


N  /  M 
+  S  (  S 
i-lX  j=l 


aiJXJ  “  bl 


where  M  may  be  any  integer  between  1  and  N. 

If  M  =  N,  we  have  the  original  problem.  For 
M  =  1,  we  have  the  simple  problem  of  minimizing 


(2.3) 


N 

+  2  (a,,x 
i«l  1 


Let  a<"> 


be  the  column  vector 


(2.4)  a(M)  - 


The  principle  of  optimality  yields  the  recurrence 
relation 

(2.5)  f*M(b)  -  min  l  X(xM  -  cM)2  +  f^_1(b  -  x^*1))] 

for  M  >  2  (cf .  [3]  ,  [4])  ,with  f2(b)  defined  by  (2.3). 

In  order  to  use  this  recurrence  relation,  we  use  the 
fact  mentioned  above  that  fM(b)  is  a  quadratic  function 
of  b, 

(2.6)  fM(b)  =  (bjQj^b)  +  2(pM,b)  +  rM, 


*111 
®2M  * 
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where  is  an  N  x  N  matrix,  pM  is  an  N-dimensional 

vector,  and  a^  is  a  scalar. 

Using  (2.5),  some  direct  algebraic  calculations 
which  we  need  not  reproduce  here  yield  the  recurrence 
relations 


(2.7) 


Si  =  St-i 

PM  =  PM-1 
rM  “  rM-l 


+ 


a(m) 

X  +  KM  ' 

(XcM  +  Pm)oc(M) 

(x  +  Km> 

(XcM  +  Pm>2 
M  ‘  +  V  ' 


where  the 

auxiliary  quantities  in  (2.7)  are  defined  by 

(2.8) 

°(M)  -  <w*(M). 

Pm  “  (Pw-la^  ^  * 

K„  =  (at»),aW), 

and 

(2.9) 

A<">  =  a<«>  ®o(M); 

here  ® 

denotes  the  Kronecker  product. 

(2.10) 

x  <5 2  y  -  ( x±y j ) i  1# j  ® 

If  we  take 

0, 


(2.11)  Qq  ®»  I,  PQ  =  0,  rQ  = 
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we  can  use  the  recurrence  relation  of  (2.7)  to  obtain 


all  Q^j,  pM,  rM  for  M  >  1.  This  saves  computing  time. 


The  minimizing  value  of  xM 

,(M) 


is 


(2.12) 


XM  = 


x+ 


m 


A  count  of  storage  requirements  and  an  estimate  of 
times  shows  that  the  foregoing  is  a  simple  feasible 
procedure  for  systems  of  dimension  up  to  90,  at  least 
with  a  computer  such  as  the  IBM-7090.  Using  various 
devices,  the  dimension  could  be  raised  considerably — at 
some  cost  in  time.  In  the  examples  that  follow,  we 
restrict  ourselves,  for  illustrative  purposes,  to  systems 
of  dimension  11. 


3.  EXAMPLE 

To  illustrate  the  method  in  practice,  consider  the 
integral  equation 

r*  1  2 

(3-1)  J  (x  —  y)^u(y)dy  =  3jj — 

0 

where  the  right— hand  side  was  chosen  so  that  the  solution 
is  u(y)  =  y. 

Let  us  use  a  quadratic  formula,  say  Simpson's  rule, 
with  11  equally  spaced  points,  y  =  0,0. 1,0. 2, . . . ,1. 
Letting  x  assume  these  same  values,  we  obtain  a  system 
of  linear  equations 


where  A  =  (a^)  and  b  =  (b1,b2, . . . are  known. 

In  Table  1,  we  show  the  results  of  various  choices 
of  c  and  X  contrasted  with  the  solution  attempted  by 
direct  methods  and  the  exact  solution. 

It  is  clear  that  some  additional  effort  is 
required  to  obtain  a  good  approximation  to  the  solution. 

Figure  1  shows  the  instability  of  the  solution 

p 

obtained  by  inversion  using  values  of  g(x)  »  x  /2  — 
2x/3  +  1/4  accurate  to  eight  significant  figures.  The 
solutions  obtained  by  means  of  dynamic  programming  are 
shown  in  Figure  2  for  various  values  of  X  using 
g(x)  accurate  to  only  three  significant  figures. 


4 .  SUCCESSIVE  APPROXIMATIONS 


To  improve  the  foregoing  results,  we  turn  first  to 
the  method  of  successive  approximations.  Consider  the 
choice  of  c  as  a  first  approximation  and  the  x 
obtained  in  this  way  as  a  second  approximation.  Generally, 
having  obtained  Xj^,  let  x^  be  determined  by  the 
condition  that  it  minimize 


(4.1)  (Ax  -  b,Ax  -  b)  +  X(x  -  Xj^qX  -  xN_1) . 

It  Is  easy  to  see  that  the  process  converges  to  a 


solution  of  Ax  =  b  for  any  initial  choice  of  c  and 
any  X  >  0,  provided  that  A  is  nonsingular. 
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Fig.  2 
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Prom  (4.1),  we  see  that 

(4.2)  (A1  A  +  M):^  -  A'b  +  XXj^. 

Since  A‘A  +  XI  is  never  singular  for  X  >  0,  we  have 

(4.3)  xN  =  (A'A  +  Xl)_1A»b  +  X(A'A  +  XI  r1^. 

Since  A  is  nonsingular,  by  assumption,  the 
characteristic  roots  of  A1 A  +  XI  are  all  greater  than 
X.  The  sequence  in  (4.3)  therefore  converges, 
geometrically,  to  a  vector  x.  From  (4.3), 

(4.4)  (A* A  +  Xl)x  »  A*b  +  XX, 

whence,  since  A‘  is  nonsingular.  Ax  °  b. 

In  the  application  of  (4.3),  the  choice  of  x  is 
crucial.  If  X  is  too  large,  the  convergence  of  (4.3) 
is  so  slow  that  computational  error  destroys  accuracy; 
if  x  is  too  small,  the  ill-conditioning  of  A»A  +  xl 
causes  grave  difficulty.  We  shall  see  how  to  overcome 
these  difficulties  to  some  extent,  below. 

5.  APPLICATION  OF  SUCCESSIVE  APPROXIMATION  PLUS  SMOOTHING 
In  Figure  3#  we  see  the  result  of  forty  iterations 
starting  with  an  initial  approximation  corresponding  to 
y/2,  with  a  value  of  X  =  0.01.  Observe  that  the 
oscillations  do  not  damp  out  in  any  strong  fashion.  If 
we  continued  the  Iteration,  round-off  error  would  soon 
eliminate  any  significance. 
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Let  us  then  smooth  the  solution  and  use  this 
smoothed  vector  as  a  new  initial  approximation.  No 
sophistication  was  used  in  the  smoothing.  Figure  4 
shows  the  great  improvement  in  accuracy  combined  with  a 
diminution  of  oscillation. 

This  process  could  be  repeated  several  times,  and 
indeed  made  part  of  the  original  program.  It  is  a 
simple  example  of  the  concept  of  sequential  computation. 

6.  EXTRAPOLATION 

Returning  to  the  problem  of  minimizing 

(6.1)  (Ax  —  b,Ax  —  b)  +  X(x  —  c,x  —  c) 

for  fixed  c,  we  know  that  the  minimizing  vector  x(X) 
approaches  xQ,  the  solution  of  Ax  =  b.  From  the 
analytic  expression  for  x(X), 

(6.2)  x(x)  =  (A'A  +  Xl)-1(b  +  Xc), 
we  see  that  as  X  -»  0, 

(6.3)  x(X)  «  xQ  +  x2X  +  XgX^  +  •••. 

Hence,  from  the  values  of  x(X)  for  small  X,  we  should 
be  able  to  estimate  xQ. 

The  point  of  the  foregoing  approach  is  that  it  is 
easier  to  compute  x(X)  for  X  >  0  because  A'A  +  XI 
is  better  conditioned  than  A.  The  basic  philosophy  is 
to  determine  the  solution  in  a  convenient  region  and  then 
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use  analytic  continuation  (here  extrapolation)  to  obtain 
the  solution  in  the  desired  region.  For  other  applica¬ 
tions  of  this  method,  see  [  5^ • 

In  Fig.  5  we  show  the  dependence  of  x(X)  for  a 
range  of  values  of  X,  where  x(X)  has  been  obtained 
using  the  foregoing  dynamic  programming  algorithm.  In 
Fig.  6,  we  show  the  relation  of  the  extrapolated 
solution,  where  again  no  sophistication  was  used  in  the 
estimation  of  xQ,  to  the  exact  solution. 

7.  COMBINATION  OF  METHODS 

Smoothing  this  approximate  solution,  we  have  an 
excellent  first  approximation  for  use  In  a  successive 
approximation  scheme.  It  is  clear  that  we  can  combine 
these  techniques  in  various  ways. 

8 .  DISCUSSION 

It  Is  generally  agreed  that  no  one  technique  will 
resolve  the  fundamental  problem  of  obtaining  sensible 
results  from  ill-conditioned  systems.  What  we  have 
wished  to  Indicate  in  the  foregoing  pages  is  that 
control  techniques  with  dynamic  programming,  successive 
approximations,  extrapolation, and  smoothing  can  yield 
worthwhile  results  In  various  cases. 

There  is  no  necessity  to  use  dynamic  programming 
algorithms  for  the  determination  of  the  minimum,  and,  in 
some  cases.  It  may  be  more  convenient  to  use  standard 


algorithms.  We  would  like  to  emphasize  that  a  dynamic 
programming  approach  has  a  built-in  stability,  as 
indicated  in  the  results  of  Sec.  3,  and  that  it  may  be 
desirable  for  this  reason  to  use  it  in  many  cases. 
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Appendix 

COMPUTING  ROUTINE 

The  computing  routine  used  in  this  Memorandum  is  given 
on  the  following  pages. 


nn 
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C  DYNAMIC  PROGRAMMING  SOLUTION 
C  MIN  (( AX-B , AX-8 ) +K < X-C , X-C >  ) 

C  QUADRATURE — SIMPSONS  RULE 

DIMENSION  CC11) *B(11>* XX (11), XY( 1 1 >  »W (11) *A< 11 t 11 ) *XI DENT (11*111 
1  XK( 1 1 )  »  X  LAM ( 12 ) »  AKRN (11.11). Q(ll, 11)  ,P(11)  * ALPHA ( 1 1  .11) 

2 . RHO ( 1 1 )  *  X  C 11)  *  G ( 1 1 ) »  XC ( 1 1  ) 

1  FORMAT  (3112) 

2  FORMAT  (E12.8) 

3  FORMAT  ( 1H1»6^X»13H MATRIX  A(I.J)///) 

4  FORMAT  ( 1 1 E 1 2 . 3  ) 

5  FORMAT  (1H0.40X.8H3  VEC TOR » 34X , 8HC  VECTOR///) 

6  FORMAT  ( 39X.E10.3.32X.E10.3  ) 

7  FORMAT  ( 1H0 » 2oX . 8HCASE  NO.  I  5 . 7X , 7HLAMBDA  =  E 1 2 . 3 / / / > 

8  FORMAT  ( 34X.2HX ( I2.3H)  =E14.5> 

9  FORMAT  (1H0.58X.17H MATRIX  ALPHA (  I  *  J  )  / / / ) 

lu  FORMAT  ( 1H0 . 39X . 10HRH0  VEC T OR » 3 3 X . 8HK  VECTOR///) 

C  NUMBER  OF  LAMBDAS  ( NLAM )  AND  MATRIX  DIMENSION  N 
READ  1.  NLAM. N,  ITER 
READ  2.  ( X LAM ( I  )  ,  1=1. NLAM) 

C  APPROXIMATE  SOLUTION  (C) 

READ  2.  (C( I ) .  1=1. N) 

DO  905  1=1, N 
XC ( I ) =C  (  I  ) 

905  CONTINUE 
C  VECTOR  B 

READ  2.  ( B ( I ) »  1=1, N) 

DO  900  1=1  , N 
G (  I  )  =8 (  I  ) 

90U  CONTINUE 

C  CALCULATION  OF  MATRIX  A 
DO  100  1=1, N 

X  I  =  I 

XX (  I  )  =  ( . l*x  I  )  -. 1 
XY ( I ) =XX( I ) 

ICC  CONTINUE 
DELTA= . L 
W(l)=  DELTA/3. 

W(N)=DELTA/3. 

M  =  N-1 

DO  101  I  =  2  »M  »  2 
W(  I )  =  (4.*DELTA)  /3 • 

101  CONTINUE 
M  =  N-2 

DO  102  1=3, M, 2 
W(  I  )  =  (2.*DELTA)/3. 

102  CONTINUE 

DO  103  1=1, N 
DO  103  J= 1 , N 

A(  I  »J) =( (XX( I )-XY( J ) )**2 )*W( J) 

103  CONTINUE 
PRINT  3 

PRINT  4.  (  ( A  (  I  , J ) »  J=1 »N)  » 1=1 »N) 

PRINT  5 

PRINT  6.  ( B ( I )  »C ( I  ) »  1=1 »  N ) 

DETERMINE  Q(l)»  P(l).  A-ND  R(l> 


nonn  n  n  n  n  nn  n  n 
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C  IDENTITY  MATRIX 
DO  104  1  =  1 » N 
DO  104  J= 1 »N 
X  I  DENT ( I » J ) =0 , 

104  CONTINUE 

DO  105  1  =  1  » N 
X IDENT (1*1  )  —  1  . 

105  CONTINUE 
C 

DO  500  I  XX  5= 1 »  NLAM 
DO  852  1  =  1  » N 
C<  I  )=XC( I ) 

852  CONTINUE 

DO  850  I  X X 7  = 1  *  I  TER 
DO  901  I  =  1  * N 
B ( I ) =G ( I ) 

901  CONTINUE 

C  KRONECKER  CROSSPRODUCT 
XK  <  1  ) =0. 

DO  106  1=1, N 

XK ( 1 ) =XK< 1 )+( A ( I , 1 ) )**2 

106  CONTINUE 
C 

DO  107  1=1, N 
DO  107  J= 1  , N 

AKRN (  I  , J ) =  ( A ( I ,1 ) *A( J ,1  )  )  / ( XLAM (  IXX5 ) +XK ( 1 ) ) 

107  CONTINUE 
C  MATRIX  Q( 1 ) 

DO  108  1=1, N 
DO  108  J= 1 » N 

G ( I , J )  =XIDENT( I »J)-AKRN( I »J) 

108  CONTINUE 

VECTOR  P( 1 ) 

DO  109  1=1, N 

P (  I  )  =  <-XLAM<  I  XX 5 )*C (  1  ) *A (  I , 1 )  ) / <  XLAM(  I XX5 ) +XK ( 1 )  ) 

109  CONTINUE 

R  <  1  I 

R  =  XLAM ( IXX5)*  <Cl 1 )**2  >-<  < XLAM (  I XX 5 ) *C ( 1 ) ) **2 ) / ( XLAM( I  XX 5 ) +XKI 1  )  ) 


DETERMINE  ALPHA ( 1 ) , RHO ( I ) 
DO  110  1=1, N 
ALPHA ( 1,1) = A (I,l| 

110  CONTINUE 

RHO ( 1 ) =0. 


BEGIN  ITERATION  OF  RECURRENCE  RELATIONSHIPS 
DETERMINE  ALPHA ( STORE  COLUMNWISE),  RHO,  AND  XK 
DO  200  IXX1=2*N 
DO  201  1=1, N 
ALPHA ( I , IXX1 ) =0. 

DO  202  J= 1 »N 

ALPHA ( I , IXX1 ) =ALPHA< I » IXX1 )+Q ( I » J  >*A( J ,IXX1 ) 
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202  CONTINUE 
201  CONTINUE 

C 

C 

RHO ( IXX1 ) *0. 

DO  203  1=1 *N 

RHO ( I XXI ) =RHO ( IXX1 )+P< I  )*A(  I, IXX1  ) 

203  CONTINUE 
C 

C 

XK ( IXX1 ) =0. 

DO  204  I = 1 »N 

XK(  I  XXI ) =  XK( I  XXI )  +  ( ALPHA (I,IXX1)*A(I,IXX1)> 

204  CONTINUE 
C 

C 

C  DETERMINE  KRONECKER  CROSSPRODUCT  OF  ALPHA 
DO  205  1=1  . N 
DO  205  J= 1 »N 

AKRN  (  I  .  J  )  =  (  (  ALPHA  (  I  .  I  XXI  )  *ALPHA( J  .  I  XXI  >  )/  (  XLAM  (  IXX5)+XK.(IXX1))) 

205  CONTINUE 
C 

C  CALCULATE  NEW  0 
C 

DO  206  1  =  1. N 
DO  206  J=1.N 
Q< I » J ) =  Q ( I . J ) -AKRN ( I .J) 

206  CONTINUE 
C 

C  CALCULATE  NEW  P 
C 

DO  207  1=1, N 

P (  I  )  =P ( I  ) —  (  (XLAM<  IXX5 )*C( IXX1 )+RHC( IXX1 )  )  *ALPHA (  I » I  XX  1  )  )  / 

1 <XLAM( IXX5 )+XK (  I XX 1  )  J 
2C7  CONTINUE 
C 

C  CALCULATE  NEW  R 
C 

R  =  R  +  XLAM ( IXX5 )*C ( IXX1 ) **2-(  ( X  LAM (  IXX5)*C(IXX1)+RH0(IXX1>  )**2 ) / 

1 (XLAM! I XX5 ) +XK ( IXX1 ) ) 

200  CONTINUE 

PRINT  7,  ( IXX5,XLAM( IXX5 J  ) 

PRINT  9 

PRINT  4  ,  ( ( ALPHA ( I » J ) »  J=1.N)»  1  =  1, N) 

PRINT  10 

PRINT  6.  ( RHO (I),XK(I)»I=1»N) 

C 

C  CALCULATE  X(N>,  X ( N-l X ( 1 ) 

C 

DO  300  IXX2=1.N 
IXX3=N+1- I XX2 

X(  I  XX  3 )  =  ( XLAM (  I XX5 ) *C (  IXX3)+RHO(  IXX3)  )/(XLAM(  IXX5)+XK(IXX3)) 

DO  301  1=1, N 

X ( I  XX 3 ) =X (  I  XX 3  )  +  (ALPHA {  I  » I XX3  >  *6 (  I  )  ) / ( XLAM (IXX5)+XK(IXX3>) 

301  CONTINUE 
C 

C  CALCULATE  NEW  VALUES  OF  8 
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DO  3C2  1  =  1, N 

B(I)=B(I)-X(IXX3)*A(I,IXX3) 
302  CONTINUE 
300  CONTINUE 

PRINT  8,  (  I » X (  I  ) ,  I =1 ,N) 

DO  851  1=1 *N 
C  (  I  )  =  X  (  I  ) 

851  CONTINUE 
850  CONTINUE 
500  CONTINUE 
CALL  EXIT 
END 


RS  X060  —  An  F402  -  Matrix  Inversion  with  Accompanying 
Solution  of  Linear  Equations. 

FAP  709/7090  Routine. 

Source:  Burton  S.  Garbow,  Applied  Mathematics  Division, 
Argonne  National  Laboratory,  Lemont,  Illinois,  dated 
February  23,  1959,  distributed  as  704  SHARE  Routine 
No .  664 . 

Revised  by  Susan  Belcher,  The  RAND  Corporation, 
September  11,  1962. 
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